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SMOOFF  -  A  FORTRAN  SMOOTHING  PROGRAM  (U) 


by 

I.H.  Grundy 


SUMMARY 

The  FORTRAN  program  SMOOFF  (SMOOthing  by  a  Filtered  Fourier 
transform)  smooths  a  "noisy"  function  of  two  variables  on  a  rectangular 
domain,  by  minimizing  an  appropriate  functional  subject  to  a  least  squares 
constraint. 
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1.  INTRODUCTION 


Techniques  for  smoothing  data  over  a  two  dimensional  region  are  well  known 
from  the  theory  of  image  processing  [1],  [2].  They  range  from  local  methods  e.g. 
using  finite  difference  schemes  such  as  five-point  Laplacian  smoothing,  to  global 
methods  in  which  a  smooth  function  is  fitted  in  some  fashion  to  the  “noisy”  input 
function  over  the  whole  region  of  interest. 

One  such  global  method  is  described  in  [2]  pp.  53-76,  in  which  the  smooth 
output  function  is  determined  by  minimising  a  given  functional,  namely  the  mean 
square  value  of  the  output  function  over  the  region,  subject  to  the  constraint  that 
the  output  function  should  be  sufficiently  close  to  the  input  function  in  a  least 
squares  sense. 

An  alternative  strategy  is  also  given  in  [2]  in  which  the  functional  to  be  min¬ 
imised  is  taken  to  be  the  mean  square  value  of  the  maximum  directional  derivative 
of  the  output  function  over  the  region.  This  article  will  describe  a  FORTRAN 
program  SMOOFF,  which  implements  this  method  on  a  rectangular  grid. 

In  SMOOFF  the  constrained  minimisation  problem  described  above  is  ap¬ 
proached  by  writing  both  the  input  function  and  the  smooth  output  function  as 
two-dimensional  Fourier  series.  The  problem  then  reduces  to  one  of  finding  the 
Fourier  coefficients  which  satisfy  the  least  squares  constraint.  This  constrained 
minimisation  problem  is  solved  using  Lagrange  multipliers. 

The  effect  of  the  minimisation  is  to  weight  the  high  frequency  coefficients  of 
the  new  Fourier  series  relative  to  the  old,  such  that  these  terms  are  attenuated. 
This  process  is  governed  by  a  single  parameter  which  enters  through  the  least 
squares  constraint.  By  varying  this  parameter  we  can  cause  as  little  or  as  much 
smoothing  as  we  desire  to  take  place. 


2.  PROBLEM  FORMULATION 

The  configuration  is  as  depicted  in  Figure  1(a).  The  grid  consists  of  N,  points 
(labelled  0,...,Nr  —  1)  evenly  spaced  a  distance  of  Aj.  apart  in  the  z-direction, 
and  Ny  points  (labelled  0, . . . ,  —  1)  evenly  spaced  a  distance  of  Aj,  apart  in 

the  y-direction.  At  each  grid  point  we  are  given  a  value  of  the  function  u(z,y). 
This  function  may  contain  a  significant  random  noise  component.  Our  task  is  to 
find  a  smooth  function  u(x,y)  which  approximates  the  noisy  input  function  u(x,y) 
sufficiently  closely. 

Later  we  will  seek  to  write  u  and  it  as  Fourier  series.  A  good  fit  of  these 
Fourier  series  to  the  actual  functions  relies,  however,  on  periodic  continuations  of 
u  and  u  outside  the  grid  being  continuous.  What  this  means  is  that  the  values  of 
u  and  it  at  opposite  edges  of  the  grid  must  be  identical,  otherwise  there  will  be 
a  step  discontinuity  in  the  periodic  continuation.  In  general  there  will  be  such  a 
discontinuity,  causing  undesirable  edge  effects  (the  Gibbs  phenomenon)  to  occur. 
This  difficulty  is  avoided  in  the  following  manner; 

We  expand  the  grid  to  form  the  one  shown  in  Figure  1(b),  The  points  in  this 
new  grid  are  labelled 

,  Vy)  =  (f  Ax,  j  Aj^),  i  =  0, . . . ,  Lx  — I,  j=0,...,Ly— 1.  (2.1) 

where  Lx  =  2yVx  -  1  and  L^  =  2/V,  —  1.  The  functions  u(x,  y)  and  u(x,y)  are  then 
extended  by  reflecting  them  in  the  lines  x  =  a/2  and  y  =  6/2,  (where  a  =  LxAx 
and  6  =  L^Aj,). 
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The  newly  constructed  functions  u  and  u  have  symmetry  in  both  the  x  and  y 
directions  and  thus  can  be  represented  in  terms  of  two-dimensional  Fourier  cosine 
series.  Further,  by  symmetry  the  values  of  u  and  u  are  equal  at  opposite  edges  of 
the  expanded  region  so  the  Fourier  representation  of  these  functions  is  excellent. 
Any  edge  effects  are  now  due  to  a  slope  discontinuity  as  opposed  to  a  function 
discontinuity,  so  their  severity  is  minimised.  (Note  that  the  last  row  and  column 
of  the  grid  are  deleted  in  Figure  1(b)  and  also  in  (2.1).  These  points  are  now 
redundant  from  the  assumed  periodicity.) 

The  criterion  for  smoothing  u  is  that  given  in  [2],  namely  that  u  should 
minimise  the  functional 


^v2 

dx^  '‘dy 


f  dxdy 


subject  to  the  constraint  that 


(2.2) 


i.-it.-i 
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1=0  y=o 


=  LiLyCr^  . 


(2.3) 


The  quantity  <t  is  the  specified  smoothing  tolerance.  If  a  is  zero,  no  smoothing 
occurs.  If  a  is  taken  to  be  too  large  the  smoothing  is  effectively  unconstrained. 
From  (2.3)  it  is  clear  that  the  optimal  value  of  a  is  closely  related  to  the  standard 
deviation  of  the  noise  on  the  data. 

The  integrand  in  equation  (2.2)  is  equal  to  the  square  of  the  maximum  of  the 
directional  derivative  at  each  point.  Thus,  a  noisy  function  with  an  underlying 
harmonic  nature  will  be  particularly  well  treated,  i.e.  without  much  distortion, 
since  the  underlying  function  minimises  I  by  definition.  This  makes  the  method 
particularly  attractive  for  dealing  with  the  output  from  SPATE,  which  has  an 
underlying  harmonic  nature  (see  (3)  and  [4]).  Even  so,  the  method  is  quite  general 
and  can  be  applied  as  a  “black-box”  smoothing  routine  for  anharmonic  functions, 
provided  the  noise  level  is  not  too  great. 

To  determine  a  suitable  u,  we  proceed  by  writing  u  and  u  as  Fourier  cosine 
series  in  x  and  y,  namely 


u 


•Yp  .Y, 

II  =rn  m=:0 


2n7rx,  2m7ry; 

- cos  — j — - 

a  b 


(2.4) 


and 


JN  p  A  , 

u  =  ^  ^  b„„,  cos 

n  =0  M»=0 


2nwx, 

- cos 
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2mwyj 
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(2.5) 


The  upper  summation  bounds  Np  and  IV,  are  chosen  so  as  to  give  a  sufficient 
number  of  terms  in  the  Fourier  expansion  to  represent  the  smoothed  function 
accurately.  Ideally  Np  =  Nr  —  I  and  IV,  =  AT,  —  1,  i.e.  the  number  of  unknowns 
is  equal  to  the  number  of  original  grid  points. 

The  Fourier  coefficients  D„,„  are  known  from  u  by  the  formulae 


Iff  I^fi  11 


4 

LsLy 
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^^u(x„yy)cos 
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ZnKXi  2mny, 

- cos - - 

a  b 
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where 


if  n  ^  0  and  m  0 
if  one  of  n  or  m  =  0 
if  n  =  m  =  0  . 


(2.7) 
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The  expression  (2.4)  for  u  is  substituted  into  (2.2)  and  (2.3).  Uniformity  of 
the  grid  is  important  here  as  it  allows  us  to  make  use  of  the  discrete  orthogonality 
property  of  sines  and  cosines,  and  this  allows  simplification  of  the  constraint  equa¬ 
tion.  The  problem  becomes  one  of  finding  the  coefficients  D„m  which  minimise 
the  sum 

N,  JV, 

f«=0  m=0 

where 

F„m  "  ^  ■52'  ’ 

subject  to  the  constraint  that 


(2.8) 


(2.9) 


JV,  N, 
n=0  m=0 


(/?„m  -  Dnmf  =  4«r^ 


(2.10) 


To  solve  this  constrained  optimisation  problem  we  use  the  Lagrange  multiplier  A 
to  write 


-V. 


AT.  JV. 


EE  ut^ n»i  (EE  )2_4ty2|=o.  (2.11) 


n=0  Ml=:0 


n=0  m=0 


Differentiation  with  respect  to  each  of  the  coefficients  in  turn  gives  that  coefficient 
in  terms  of  its  known  hatted  counterpart,  i.e. 


^nm  — 


ADnrM 

A  +  F„, 


(2.12) 


for  all  n  and  m. 

The  next  step  is  to  find  the  single  unknown  A.  This  is  accomplished  by 
substituting  (2.12)  into  the  constraint  equation  (2.10),  to  give  a  single  non-linear 
algebraic  equation  for  A,  namely 


A-,  .Y, 


ip2 

t*w^nin  f\2 
2  w 


^MW*) 

ri=0  wi=0  '  ' 


4(t^  =  0  . 


(2.13) 


This  equation  is  solved  numerically  by  a  Newton  iteration  scheme.  We  observe 
that  £■  is  a  monotonically  decreasing  function  of  A,  hence  (2.13)  will  have  a  unique 
solution  for  A.  A  non-negative  solution  is  guaranteed  provided  £(0)  >  0.  The 
critical  case  occurs  when  £{0)  =  0.  This  case  corresponds  to  all  the  Fourier 
coefficients  of  u  being  zero  (see  2.12),  with  the  exception  of  the  constant  term 
Don-  That  is,  it  corresponds  to  unconstrained  smoothing  and  total  loss  of  the 
original  information.  Using  (2.13)  we  can  determine  the  critical  a  value,  namely 
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yp  A*. 


(2.14) 


which  gives  rise  to  this  unconstrained  smoothing. 

With  A  finally  determined,  it  remains  to  calculate  the  Fourier  coefficients  D„,„ 
from  equation  (2.12),  and  to  substitute  these  back  into  the  equation  (2.4)  for  u. 
From  this  series  u  can  be  evaluated  anywhere  within  the  rectangular  region,  and 
in  particular  at  the  original  grid  points. 


3.  PROGRAM  DESCRIPTION 

The  smoothing  algorithm  described  above  has  been  implemented  in  the  double 
precision  FORTRAN  subroutine  SMOOFF.  This  subroutine  is  divided  into  three 
distinct  stages. 

The  first  stage  involves  expanding  the  grid  and  determining  the  Fourier  co¬ 
efficients  DH(N,M)  of  the  expanded  UDATA  on  this  grid.  Symmetry  of  UDATA 
about  the  lines  x  =  a/2  and  y  =  i/2  is  used  to  offset  the  extra  computational 
effort  caused  by  the  grid  expansion. 

The  second  stage  of  SMOOFF  involves  determining  the  Lagrange  multiplier 
ALAM  by  Newton  iteration.  This  is  the  fastest  stage,  being  an  order  of  magni¬ 
tude  quicker  than  stages  one  and  three.  Convergence  always  occurs  within  ten 
iterations.  The  critical  SIG  value  SIGMAX  is  also  determined  here.  If  SIG  is 
greater  than  SIGMAX  then  the  subroutine  exits  immediately  without  going  on  to 
the  third  stage. 

The  third  stage  of  SMOOFF  involves  back-substitution  of  ALAM  to  give 
the  Fourier  coefficients  DS(N,M)  of  the  output  function  UOUT.  Further  back- 
substitution  gives  the  smooth  output  UOUT  itself  at  the  grid  points.  Again  sym¬ 
metry  is  used  to  reduce  the  computational  effort  by  something  approaching  a  factor 
of  four. 

A  typical  call  to  SMOOFF  looks  like 

CALL  SMOOFF (NX . NY . DX . DY . SIG .SIGMAX . NA . UDATA . UOUT . lOP . lER . IMESSG) . 
The  input  parameters  are 

NX  the  number  of  original  grid  points  in  the  i-direction, 

NY  the  number  of  original  grid  points  in  the  y-direction, 

DX  the  grid-spacing  in  the  s-direction, 

DY  the  grid-spacing  in  the  y-direction, 

SIG  the  smoothing  tolerance, 

NA  the  size  of  UDATA  and  UOUT  in  the  calling  program, 

UDATA  the  array  of  raw  function  values,  and 
lOP  the  fast/slow  run  option 

The  input  parameter  lOP  deserves  a  special  mention  here.  This  parameter  controls 
the  number  of  Fourier  terms  and  JV,  which  are  considered.  The  default  is  for 
the  ideal  case:  Nf,  =  —  i  and  =  Ny  —  I.  For  lOP  equal  to  unity  however, 

N/,  and  N,j  are  halved  from  their  ideal  values,  in  effect  setting  the  remaining  (high 
frequency)  tenns  to  zero  and  truncating  the  Fourier  series.  This  is  done  in  order 
to  reduce  the  computation  time  by  a  factor  of  four.  A  side  effect  though,  is  that 
the  smoothing  parameter  cr  no  longer  relates  directly  to  the  noise  level  as  in  (2.3). 
A  significantly  smaller  value  of  <t  needs  to  be  chosen,  a  fact  which  is  clear  from 
inspection  of  (2.11)  and  (2.12).  Nevertheless,  for  cases  with  a  lot  of  grid  points, 
say  greater  than  40  x  40  the  fast  run  option  gives  reasonable  results  in  acceptable 
time.  The  option  should  not  be  used  however,  if  a  derivative  of  u  is  required,  or  if 
u  itself  is  expected  to  contain  high  frequency  components. 

The  output  parameters  of  SMOOFF  are 

SIGMAX  the  critical  value  of  SIG, 

UOUT  the  array  of  smoothed  function  values, 

lER  the  error  number,  and 

IMESSG  the  error  /  nomial-completion  message. 
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To  use  SMOOFF  the  following  declarations  should  be  made  in  the  calling  program: 

INTEGER  NX.NY.IER.IOP 

CHARACTER*50  IMESSG 

REAL*8  DX.DY.SIG.SIGHAX 

DIMENSION  UDATA (0 : NA, 0 : NA) .UQUT(0 : NA , 0 : NA) 

PARAMETER  (NA>199} . 


4.  SAMPLE  RESULTS 


Sample  results,  consisting  of  surface  and  contour  plots  obtained  using  the 
NCAR  plotting  package,  are  presented  here  for  the  simple  test  function  u(x,j/)  = 
zy,  both  with  and  without  added  random  (Gaussian)  noise.  The  results  are  pre¬ 
sented  for  a  uniform  50  x  50  grid  on  |0, 1]  x  [0,  Ij.  In  each  case  the  calculation  of  the 
smoothed  approximation  took  around  110  CPU  secs  using  the  slow  run  option,  of 
which  about  half  was  spent  in  stage  1  and  about  half  in  stage  3.  (For  a  120  x  120 
grid  the  total  CPU  time  is  of  the  order  of  2  hrs.) 

Figure  2(a)  shows  a  surface  plot  of  the  original  function  u  without  noise. 
Figure  2(b)  shows  the  corresponding  contour  plot.  The  contours  are  (naturally) 
the  hyperbolae  zy  =  c. 

Figure  3(a)  shows  ti  with  added  Gaussian  noise  with  a  standard  deviation  of 
0.05.  Figure  3(b)  shows  the  corresponding  contour  plot.  The  original  contours  are 
now  largely  obscured  by  noise. 

Figure  4(a)  shows  the  effect  of  SMOOFF  on  the  noisy  u  of  Figure  3.  The 
smoothing  level  is  =  0.05  (ffmax  ~  0.222),  which  is  precisely  equal  to  the  noise 
level.  Clearly,  this  is  very  near  the  optimal  value  for  a. 

Finally,  Figures  5(a)  and  5(b)  show  the  effect  of  an  excessively  high  choice 
of  the  smoothing  tolerance  (er  =  0.2).  There  is  significant  loss  of  information, 
demonstrated  by  gross  distortion  of  the  contours  and  flattening  of  the  surface.  (It 
should  be  pointed  out  here  that  the  surface  plot  is  not  drawn  to  scale,  and  is  rather 
flatter  than  it  appears!) 


5.  CONCLUSION 


The  FORTRAN  subroutine  SMOOFF  implements  the  method  of  smoothing 
noisy  data  on  a  rectangular  grid,  described  in  [2).  Numerical  evidence  presented  for 
sample  data  in  Section  4  shows  that  SMOOFF  is  very  effective  provided  sufficient 
care  is  taken  with  the  size  of  the  smoothing  parameter.  Applications  have  already 
been  found  at  ARL  in  smoothing  output  from  SP.VTE  (see  [3],  [4]). 
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APPENDIX  -  PROGRAM  LISTING 


C 

ccccccccccaxxxxxxxxxxccaxcccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTI NE  SMOOFF (NX . NY . DX . OY . SIG . S IGNAX . NA . 

*  UDATA.UOUT.IOP.IER.IMESSG) 

C 

cccccccccccccccccccccccccocccc<xxx;cccccccccccccccccc<x;ccccccccccccc 

AUTHOR:  Ian  Grundj  12/2/87. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


This  subroutine  snooths  noisy  input  data  UDATA,  defined  on 
a  rectangular  region  containing  an  NX*NY  distribution  of 
points  a  distance  of  DX  and  DY  apart  respectively.  (There 
is  no  restriction  on  the  oddness  or  evenness  of  NX  and  NY.) 
The  resulting  smoothed  data  are  written  into  UOUT.  The 
degree  of  smoothing  is  controlled  by  the  parameter  SIG. 

This  process  is  implemented  in  three  stages.  In  Stage  1, 
the  grid  and  the  raw  input  data  UDATA  are  reflected  in 
vertical  and  horizontal  lines  of  symmetry  to  form  a  grid  of 
size  four  times  the  original.  UDATA  is  then  expressed  as  a 
two-dimensional  Fourier  cosine  series  whose  coefficients  DH 
are  then  computed.  Here  the  number  of  grid  points  in  the  x 
direction  is  LX=2*NX-1  and  the  number  in  the  y  direction  is 
LY*2*NY-1  (since  periodicity  is  assumed  at  the  far 
boundary) . 

The  Fourier  coefficients  are  input  to  Stage  2,  where  the 
Lagrange  multiplier  ALAM  is  used  to  satisfy  a  variational 
principle  (minimise  the  integral  of  the  squares  of  the  first 
derivatives  of  UOUT  over  the  region),  subject  to  the 
constraint  that  the  "distance”  between  UOUT  and  the  known 
Fourier  series  for  UDATA  is  smaller  the  smoothing  tolerance 
LX*LY*SIG*SIG. 

ALAN  is  input  to  Stage  3  where  the  Fourier  coefficients  DS 
of  UOUT,  and  finally  UOUT  are  computed. 

NOTES: 

(1)  SIG=0  implies  that  UOUT  is  exactly  equal  to  the 
Fourier  series  for  UDATA,  l.e.  effective  smoothing 
will  not  take  place. 

(2)  If  SIG  is  taken  too  large,  i.e.  larger  than  SIGMAX 
which  is  computed  within  STAGE  2,  the  smoothing  is 
effectively  unconstrained,  ALAN  is  set  to  zero  and  the 
smoothed  function  UOUT  becomes  constant  over  the  region. 

(3)  Error  conditions. 

IER=0  No  Error 

IER=1  Fatal  Error 

IER=2  Exceptional  Cases  (i)  and  (2) 

INESSG  gives  the  error  message. 
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c 

C  (4)  I0P=1.  Time  saver  mode.  The  program  runs  for  a  reduced 
C  number  of  Fourier  coefficients.  In  effect  setting  the 

C  remainder  to  zero. 

C 

C  DECLARATIONS : 

C  In  the  calling  program  the  variables  should  be  declared  in 

C  the  following  manner ; 

C 

C  INTEGER  IIX.NY.IER.IOP 

C  CHARACTER*50  IMESSG 

C  REAL*8  DX.DY.SIG.SIGMAX 

C  DIMENSION  UDATA(0:NA.0:NA) .U0UT(0:NA.0:NA) 

C  PARAMETER  (NA=199) 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

INTEGER  LX , LY . NX . NY . NA . lER . lOP 
CHARACTER*50  IMESSG 
REAL » 8  DX.DY.SIG.SIGMAX.ALAM 
DIMENSION  DH(0;NU.0:NU) ,DS(0;NU.0:NU) . 

*  UDATA(0:NU.O.-NU).UOUT(0:NU.O;NU)  , 

*  DSNX(0:NU2) .DCSX(0:NU2) .DSNY(0:NU2) .DCSY (0 : NU2) 
PARAMETER  (NU=199 . NU2=399) 

C 

PI=4 . 0D0»DATAN ( 1 . ODO) 

LX=2*NX-1 

LY=2*NY-1 

A=LX*DX 

B=LY*DY 

TNN=16 . ODO/DFLOAT (LX*LY) 

PNX=2.0D0*PI/DFL0AT(LX) 

PNY=2.0D0*P1/DFL0AT(LY) 

C 

C  If  I0P=1  the  program  works  with  a  reduced  number  of  Fourier 

C  coefficients,  in  effect  setting  the  remainder  to  zero.  This 

C  becomes  attractive  beyond  NX=40,NY=40. 

C 

IF  (lOP.EQ.l)  THEN 
NXH=(NX-l)/2 
NyH=(Ny-l)/2 
ELSE 

NXH=NX-1 

NyH=NY-l 

ENDIF 

C 

C  Error  Tests.  If  there  is  an  error,  IER=1. 

C 

IER=0 

IF  (DX.LE.O.OR.DY.LE.O)  THEN 

IMESSG=’DX  or  DY  is  negative  or  zero’ 

IER=1 

ENDIF 

IF  (NX.LE.O.OR.NY.LE.O)  THEN 

IMESSG='NX  or  NY  is  negative  or  zero’ 

IER=1 

ENDIF 

IF  (SIG.LT.O.ODO)  THEN 
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oo  ooo  oooQoonooooo  ooo 


t  IMESSG= 'Smoothing  tolerance  SIG  must  be  non-negative' 

'  IER=1 

vENDIF 

aF  (NA.HE.NU)  THEN 

'  IHESSG='Array  dimensions  NA  and  NU  must  be  equal’ 


1  IER=1 
HiNDIF 

(lER.EQ.l)  RETURN 
Tipt  for  zero  sigma 


IF  \(SIG.LE.1E-10)  THEN 
DO  2  1=0, NX-1 
V,D0  1  J=0.HY-1 
'\  UOUT(I.J)=UDATAa.J) 

1  CONTINUE 

2  CONTINUE 

IHKSG='SIG  is  zero.  No  smoothing  done.' 

IE^2 
RETBRN 
ENDIF  I 

:cccccccccc^ccccccccccccccccccccccccx:cccccccccccccccccccccccccccc 

STAGE  1\: 

Takes  rdjw  data  UDATA  and  computes  its  two-dimensional 
Fourier  aeries  with  coefficients  DH(N.M). 

xcccccccccccAccccccccccccccccccccccccccccccccccccccccccccccccccc 

Set  up  ar|ay  of  sines  and  cosines 

\  DO  101  I=0|LX 

DSNX(I)=-»IN(PNX*I) 

DCSX(n»M0S(PNX*I) 

101  CONTINUE  \ 

DO  102  >0.1X 

DSNy(J)-DSlN(PNY»J) 

DCSY(J)-DC«(PNY»J) 

102  CONTINUE  I 

Find  DH  for  am  N  and  N. 

DO  110  N-0.NXH1 
DO  109  H-O.NW 
DH(N,M}-UDiCTA(0.0) 

DO  106  I-l.lx-1 
NI-MOD(N*iLX) 

DH(N.M)-DHIN.M)+2.0D0*DCSX(NI)*UDATA(I.0) 

106  CONTINUE  \ 

DO  106  J-l.NTfrl 
MJ=MOD(M*J,Iy) 

DH (N . M) =DH (1 . M) +2 . 0D0*DCSY (M J) »UDATA (0 , J) 

106  CONTINUE  \ 

DH(N.M)-DH(N,m1/4.0D0 

Sum  over  space  varnbles 


DO  108  I=1.MX-1 
TD=O.ODO 
DO  107  J=1.NY-1 
MJ=MOD(M*J.LY) 

TD=TD+UDATA(I . J) *DCSY (M J) 

107  CONTINUE 
NI=MQD(N»1.LX) 

DH (N , M) =DH (N . M) +DCSX (N I) *TD 

108  CONTINUE 
C 

WT=2.0D0 

IF  (N.GT.O.AND.M.GT.O)  WT=1.0D0 
IF  (N.EQ.O.AND.H.EQ.O)  WT=4.0D0 
DH (N . M) “DH (M . M) *TNM/WT 

109  CONTINUE 

110  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  STAGE  2: 

C 

C  Finds  Lagrange  multiplier  ALAN  on  which  the  smoothed 
C  solution  depends.  This  involves  solving  a  non-linear 
C  algebraic  equation  for  ALAN  using  a  modified  Newton 
C  iteration  scheme. 

C 

CCCCCCCCCCCCOCCCCCCCCCCCCOCOXICCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  Find  initial  estimate  for  ALAN  and  critical  sigma  value 
C  SIGNAX. 

C 

ALAN“0.0D0 
SIGNAX-O.ODO 
DO  202  N“0.NXH 
DO  201  N>0.NYH 

IF  (N.EQ.O.AND.N.EQ.O)  GOTO  201 

FNN=N*N/A/A+N*N/B/B 

WT“2.0D0 

IF  (N.GT.O.AND.N.GT.O)  WT=1.0D0 
SNN“WT*DH(N,N)*DH(N.N) 

SIGNAX“SIGNAX^SNN 

HNN>SNN*FNN*FNN 

ALAN-ALAN+HNN 

201  CONTINUE 

202  CONTINUE 

SIGNAX-0 . 5D0*DSqRT (SIGNAX) 

C 

C  Test  for  unconstrained  case 
C 

IF  (SIG.GE. SIGNAX)  THEN 
DO  204  I-O.NX-1 
DO  203  J“0,NY-1 
U0UT(I.J)“DH(0.0) 

203  CONTINUE 

204  CONTINUE 

INESSG“*SIG  is  too  large.  Smoothing  unconstrained.' 

IER=2 

RETURN 

ENDIF 

C 
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C  Iterate  to  find  ALAM 
C 

IF  (SIG.LT.0.2SD0*SIGMAX)  THEN 
ALAM=0 . 6D0* (DSQRTCALAM) /2 . ODO/SIG) 

ELSE 

ALAM=O.ODO 

ENDIF 

205  CONTINUE 
FLAN-O.ODO 
FDLAM-O.ODO 

DO  207  N=°0,NXH 
DO  206  M^O.NYH 

IF  (N.EQ.O.AND.N.EQ.O)  GOTO  206 

FNM»N*N/A/A+K*M/B/B 

WT-2.0D0 

IF  (N.GT.O.AND.M.GT.O)  Vrr=1.0D0 
HNM=WT»FNM*FMM*DH(H.M)*DH(N.M) 

AF=(ALAM+FNI0 

FLAM-FLAM+HNM/AF/AF 

FDLAM-FDLAM+HMM/AF/AF/AF 

206  CONTINUE 

207  CONTINUE 
FDLAM=-2.0D0*FDLAM 
FLAM=FUM-4 . 0D0*SIG»SIG 
ALAMOLD»ALAM 
ALAM-ALAM-FLAM/FDLAM 

C 

C  If  ALAN  becomes  negative,  halve  the  initial  guess  for  it. 

C 

IF  (ALAM.LT.O.ODO)  ALAN  »  0.SD0*ALAM0LD 
C 

IF  (ABS<(ALANOLD-ALAM)/ALAN) .GT.IE-IO)  GOTO  20S 
C 

ccccccccccccoccccccccccccocaxccccccccccccccccccccccccccccccccccccc 

c 

C  STAGE  3; 

C 

C  Substitutes  the  calculated  value  for  ALAN  into  the 
C  expressions  linking  the  Fourier  coefficients  DH  and  DS  for 
C  all  N  and  N.  UOUT  is  then  calculated  at  every  grid  point 

C  using  its  now  known  Fourier  series. 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  Calculate  the  smoothed  Fourier  coefficients  DS. 

C 

DO  302  N=0,HXH 
DO  301  N-0,NYH 

IF  (N.EQ.O.AND.N.EQ.O)  GOTO  301 

FNN-N*N/A/A+N*N/B/B 

AF-ALAN+FNM 

DS (N . M) =ALAH*DH (N . M) /AF 

301  CONTINUE 

302  CONTINUE 
DS(0.0)=DH(0.0) 

C 

C  Use  the  calculated  coefficients  to  calculate  UOUT  at  all 
C  grid  points. 

C 

DO  306  1=0. NX- 1 
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DO  305  J=0.MY-1 
C 

UOUT(I,J)»O.ODO 
DO  304  N-O.NXH 
TD-0 . ODO 
DO  303  N-O.NYH 
M>NOD(N*J,LY) 

TD-TD+DS (N . M) *DCSY (M J) 

303  CONTINUE 
NI-MOD(N*I.LX) 

UOUTd .  J)»UOUT(I .  J)+DCSX(NI)  *TD 

304  CONTINUE 

305  CONTINUE 

306  CONTINUE 
C 

IMESSG> ' SNOOFF  Normal  Completion.  Have  a  nice  Day.' 

RETURN 

END 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCOCCCCCCCCCCCOCCCCCCCCCCCCOCCCCCCCCCCCCCCCC 


Original  grid. 


Expanded  grid. 


FIGURE  2(a).  Surface  plot  of  original  input  function. 


FIGURE  3(a).  Surface  plot  of  noisy  input  function  (original  input  with  added 
Gaussian  noise  -  standard  deviation  O.OSj. 


FIGURE  4(a).  Surface  plot  of  noisy  input  function  after  application  of  SMOOFF 

at  a  =  0.05. 
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FIGURE  S(a).  Surface  plot  of  noisy  function  after  processing  by  SMOOFF 
showing  the  effect  of  an  excessive  ex  value  (a  =  0.2). 
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